Two-dimensional (2D) Dirac materials,^[1] characterized by linear energy dispersion near the Fermi level, have shown potential applications in next generation nanoscale devices,^[2,3] due to their unique electronic properties such as ballistic charge transport,^[4] high carrier mobility,^[5] and quantum Hall effect.^[6] Due to the rigorous requirements for the symmetry, parameters, Fermi level, and band overlap in the materials to display Dirac cone features,^[7] 2D Dirac materials are relatively rare and limited to very few types of lattices, among which the honeycomb lattice provides the most common platform. Graphene, as an archetypal 2D Dirac material with honeycomb lattice, has been intensively studied during the past decade. In graphene, all the carbon atoms lie in the same plane, and the active bands near the Fermi level arise from π bonding and are composed of the p_z orbitals directly normal to the graphene plane. The other two in-plane p orbitals (p_x,y) exhibit both orbital degeneracy and spatial anisotropy and they hybridize with the s bands, forming fully filled sigma-bonding bands. Similar views work for graphene analogues including silicene, gemanene, and stanene,^[8–13] except that the honeycomb lattice is buckled and all the atoms do not lie in the same plane any more. But still, the Dirac states in the free standing structures are composed of the p_z orbitals. However, in practice, free standing forms of graphene analogues have not yet been made so far. On the other hand, the relatively strong coupling between p_z orbitals and the substrates makes the p_z orbital Dirac states in graphene analogues unobservable.^[14,15] Alternatively, recent research has proposed to realize p_x,y-orbital Dirac states in silicene grown on SiC substrate.^[16] Due to the in-plane nature of p_x,y orbitals, the coupling from the substrate is expected to have less effect on this kind of Dirac states, and thus p_x,y-orbital Dirac states could be more robust for practical applications. In fact, just as the p_z orbital, p_x,y orbitals often play an important role in the honeycomb lattice and they have been systematically investigated in the contest of ultracold-atom optical lattices.^[17–23] Recently, Dirac states composed of p_x,y orbitals have been reported in several 2D systems including fluoridated tin film,^[10,24,25] functionalized germanene,^[26] functionalized group V monolayers BiX/SbX (X = H, F, Cl, Br),^[27,28] and some organic materials.^[29–31] Besides, a variety of topological band structures can be realized by p_x,y orbitals in the honeycomb lattice.^[27,32] Owing to the emerging interest and potential importance, the p_x,y-orbital Dirac states in the honeycomb lattice are required to be comprehensively understood, i.e., to obtain better theoretical insight into them through low-energy-band models.

To this end, we first construct a four-band tight-binding model for the p_x,y-orbital Dirac states considering not only the nearest neighbor hopping interactions but also the buckling effect of the atomic planes, as the honeycomb lattice composed of p_x,y-orbital Dirac states is often buckled. Unlike p_z-orbital Dirac states, which contain only two linearly dispersed bands near the Dirac cones, we find that p_x,y-orbital Dirac states also contain two additional narrow bands that are flat in the limit of vanishing π bonding. The results are in agreement with previous studies.^[17,19,32] Most importantly, we analytically obtain the linear dispersion relationship between the energy and momentum vector near the Dirac cone. We find that the Fermi velocity is determined by the hopping interactions through both the π bonding and the σ bonding of p_x,y orbitals, which is in contrast to the case of p_z-orbital Dirac states. Consequently, p_x,y-orbital Dirac states show more flexibility for engineering, since the changes of Fermi velocity are more sensitive to the changes of lattice constant and buckling angle. We further confirm the validation of our tight-binding model by direct first-principles computations of model material systems including hydrogenated monolayer Bi and Sb honeycomb lattices. Our work conduces to the understanding of p_x,y-orbital Dirac states in honeycomb lattices, which will also be useful for the application of this family of materials to nanoelectronics.

We perform the first-principles calculations using the density functional theory (DFT),^[33,34] as implemented in the VASP code.^[35–38] The interaction between electron and core is included using the frozen-core projected augmented wave (PAW)^[39,40] approach, and the generalized gradient approximation (GGA) formulated by Perdew, Burke, and Errnzerhof (PBE)^[41,42] is adopted. The structures are relaxed until the residual force on each atom is less than 0.01 eV/Å and total energy converges to 10⁻⁸ eV, with the cutoff energy for plane-wave basis functions set to be 400 eV. By neglecting the second and higher order term with respect to q², the Fermi velocity is estimated for BiH and SbH by fitting the Dirac bands at k = K + q to the expression v_F = E(q)/ħ|q|. For the discussion of the semimetallic Dirac states, the spin–orbital coupling effect is not considered in this work as it will break the Dirac states and open bandgaps.

As shown in Figs. 1(a) and 1(b), the graphene-like honeycomb lattice structure can be seen as two triangular sublattices A and B. It effectively has two atomic sites per unit cell and the lattice vectors can be written as 1 where is the length of the lattice vector. For a honeycomb lattice in the absence of buckling, the atoms of sublattices A and B lie in the same plane; however, for the buckled case, the atoms in the sublattice A lie in a parallel atomic plane with those in the sublattice B (see Fig. 1(a). To quantify to what extent the honeycomb lattice is buckled, here we define the angle between the atomic bonds and the unbuckled atomic plane as buckling angle θ, as seen in Fig. 1(c). For p_x,y-orbital Dirac states, each atomic site has two orbitals, p_x and p_y. Thus electrons can hop from one orbital at one atomic site to two orbitals at a nearest atomic site. By just considering the nearest hopping terms, the tight-binding Hamiltonian for p_x,y electrons in a graphene-like honeycomb structure is 2 where a_i,μ annihilates (creates) an electron with orbital μ (μ = p_x, p_y) at site R_i in the sublattice A (an equivalent definition is used for sublattice B), and t_i,j,μ,ν are the nearest-neighbor hopping integrals. We choose as a basis set as depicted in Fig. 1(a), to construct the Hamiltonian matrix. The resulting Hamiltonian in the k-representation is 3 where x and y denote p_x and p_y orbitals, while A and B denote sublattice sites, respectively. The matrix elements are 4 where are the hoping integrals between different orbitals at different atomic sites, connected by δ_i, as indicated in Fig. 1(a); δ₁, δ₂, and δ₃ are the three nearest-neighbor vectors in real space and given by 5 According to the Slater–Koster method,^[43] the hoping integrals have the following relationships: 6a 6b 6c 6d 6e 6f Note that V_ppπ and V_ppσ are the hopping integrals through the in-plane π and σ bonds, respectively. In the buckling case, the hopping integrals will become the superposition of in-plane π and σ bond integrals, i.e., cV_ppσ+(1−c)V_ppπ, where c is equal to cos²θ.

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) (a) Top and side view of the buckled honeycomb lattice, in which hopping integrals between atomic orbitals are indicated with . Sublattice A (B) is indicated with red (blue) balls. The lattice vectors a1, a2, and the nearest distances between sublattices δ1, δ2, δ3 are also marked. (b) Sublattice A with its nearest neighbours B’s. (c) Diagram showing definition of buckling angle θ.

Considering all the above equations, the solved eigenvalues as functions of momentum vector for the tight-binding Hamiltonian of p_x,y orbitals in honeycomb lattices are 7a 7b where the plus and minus signs refer to the upper and lower bands, respectively. Note that unlike the case of p_z-orbital Dirac states in graphene where the eigenvalues only depend on the hopping integrals through π bonds, the eigenvalues for the p_x,y-orbital Dirac states in the honeycomb lattice depend on both π bond hopping integrals (V_ppπ, which is positive according to Fig. 1(a)) and σ bond hopping integrals (V_ppσ, which is negative according to Fig. 1(a)).

To obtain a better understanding of the dispersion relationships between energies and momentum vectors, in Figs. 2(a)–2(c), we draw the band structures in the whole reciprocal space for the no-buckling case with different choices of V_ppπ and V_ppσ. We chose a rectangular region, which includes the first Brillion zone (FBZ), to show the 3D energy dispersion. As clearly seen, Dirac cones emerge in all cases. Besides, the upper and lower bands are symmetric with respect to the Fermi level. However, such Dirac states in real systems are usually not purely composed of p_x,y orbitals but are hybridized with other orbitals, therefore, the four bands are not rigorously symmetric. For example, we calculate the orbital-projected band structure of hydrogenated monolayer Bi and Sb honeycomb lattices at their equilibrium lattice constants in Figs. 3(a)–3(d). The equilibrium lattice constants and buckling heights for hydrogenated monolayer Bi and Sb honeycomb lattices in the calculations are 5.53 Å, 0.08 Å, and 5.29 Å, 0.08 Å, respectively. We can see that while the Dirac states are indeed mainly composed of p_x,y-orbital bands and keep the symmetric feature, they also overlap with other bands, making the Dirac bands not rigorously symmetric any more. However, if sufficiently large strain is applied to hydrogenated monolayer Bi and Sb honeycomb structures, to ensure that the p_x,y-orbital bands are isolated from other bands, we do observe four clearly symmetric bands around the Dirac states, which are shown in Fig. 4. The magnitude of strain is described by (a − a₀)/a₀, where a₀ and a denote the lattice constants of the unstrained and strained systems, respectively.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Band structures of px,y-orbital Dirac states for the unbuckled case with different choices of hopping integrals of Vppπ and Vppσ: (a) Vppπ = 0.33 eV and Vppσ = −0.97 eV; (b) Vppπ = 1.30 eV and Vppσ = 0.00 eV; (c) Vppπ = 0.00 eV and Vppσ = −1.30 eV. Note that, the flat bands appear when Vppπ = 0 or Vppσ = 0.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) Bi (a) py and (b) px orbital projected electronic band structure for the BiH buckled monolayers from first-principles calculations. Sb (c) py and (d) px orbital projected electronic band structure for the SbH buckled monolayers from first-principles calculations. Fermi energy is set to be zero. Symbol (colored dots) size is proportional to the contribution of the corresponding state. Red dots are for the Bi or Sb py orbitals and green dots for the Bi or Sb px orbitals. Note that the four bands near the Fermi level are mostly from the px and py orbitals.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Band structure of (a) BiH and (b) SbH under different strains Δa/a: 20%, 40%, 60%. The px,y-orbital bands are isolated from s and pz bands. The symmetric feature of bands with respect to fermi level becomes clearer, when strains increase. A four-orbital model becomes more accurate for the structure under big strains.

Another interesting feature is that when V_ppπ = 0 (or V_ppσ = 0), two flat bands corresponding to (or ) emerge at the upper and the lower parts of the Dirac cones, as shown in Figs. 2(b) and 2(c). With charge carriers of high effective mass and low mobility, these two flat bands could offer interesting features^[44] such as spatially confined electrons, high peaks in the density of states (DOS), and Lifshitz transitions^[45] (if the two bands are located at the Fermi level). In a practical system, such a band nature can be realized in honeycomb lattices with very large effective lattice constants, i.e., in metal-organic-frameworks, in which the π bonding is much weaker than the σ bonding.^[31] We also demonstrate this feature by using our model systems through stretching the lattices. As shown in Fig. 4, the two bands at the upper and the lower parts of the Dirac cones become flatter and flatter with tensile strain increasing. Our findings are in agreement with previous results.^[17,19,32]

Now we turn our attention to discuss the properties of the Dirac cones. Near the corner points of the Brillouin zone (BZ), i.e., near at k = K + q with |q| ≪ |K|, we obtain the linear dispersion relationship between energy and momentum by expanding the full band structure, Eq. (7a), i.e., , where q is the momentum measured relatively to the Dirac points and v_F is the Fermi velocity which is given by or . Note that the Fermi velocity adopts similar dependence on the hopping parameter and lattice constant in the case of the p_z-orbital Dirac state, where the Fermi velocity for the unbuckled honeycomb structures is v_F = 3/2 at with t being the hopping parameter between nearest p_z orbitals. The different thing is that in the case of p_x,y-orbital Dirac states, Fermi velocity is determined by the summed hopping strength of σ bonding and π bonding while only π bonding hopping strength plays a role in the case of p_z-orbital Dirac states. The consequence is that in the case of p_x,y-orbital Dirac states, the changes of Fermi velocity will be much more sensitive to the changes of lattice constant and buckling angle, which could offer an advantage for flexible electronics and mechanical sensors. Indeed, as shown in Fig. 5(a), we find that the relative Fermi velocity changes in BiH and SbH are larger than in graphene when the lattice constants are stretched by the same ratio. When the lattice constants are compressed, the relative Fermi velocity change in BiH is still larger than in graphene. However, the relative Fermi velocity change in SbH turns smaller with larger compressive lattice constant ratio. This is because more and more p_z orbitals are involved in the Dirac states of SbH when the lattice constants are compressed, as shown in Table 1. Figure 5(b) shows the plots of relative Fermi velocity changing with buckling angle for graphene, BiH, and SbH. To ignore the changes of the σ and π bond integrals, we keep the bond length fixed while changing the buckling angle θ. Because v_F ∼ acos²θ and a ∼ cosθ, v_F is linearly dependent on cos³θ, as seen in Fig. 5(b). We can see that the relative Fermi velocity changes in both buckled BiH and SbH are larger than in buckled graphene when the buckling angles are changed by the same amount, demonstrating that p_x,y-orbital Dirac states do show more engineering sensitivities and flexibilities than p_z-orbital Dirac states.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Changes of relative Fermi velocity variation in a unit of the Fermi velocity (v0) without strain in BiH, SbH, and graphene with change of (a) strain Δa/a0, and (b) cos3θ.

Table 1.

Table 1.

Table 1. Orbital contribution of Sb pz orbital in SbH layer at each strain. More pz orbitals are involved in the Dirac states when the lattice constants are compressed. Difference in contribution of pz orbital makes the change of relative Fermi velocity smaller with larger compressive lattice constant ratio. . Δa/a0 −0.06 −0.05 −0.04 −0.03 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06 pz 0.027 0.025 0.023 0.021 0.020 0.019 0.018 0.017 0.016 0.015 0.014 0.014 0.013

Table 1. Orbital contribution of Sb pz orbital in SbH layer at each strain. More pz orbitals are involved in the Dirac states when the lattice constants are compressed. Difference in contribution of pz orbital makes the change of relative Fermi velocity smaller with larger compressive lattice constant ratio. .

In this work, we have constructed a four-band tight-binding model for the p_x,y-orbital Dirac states considering not only the nearest neighbor hopping interactions, but also the buckling effect. We find that p_x,y-orbital Dirac states contain two additional narrow bands, which are flat in the limit of vanishing π bonding. The results are in agreement with previous studies. Most importantly, we analytically obtain the linear dispersion relationship between energy and momentum vector near the Dirac cone. We find that the Fermi velocity is determined not only by the hopping through π bonding but also by hopping through the σ bonding of p_x,y orbitals, which is in contrast to the case of more common p_z-orbital Dirac states. As a result, p_x,y-orbital Dirac states show more flexible engineering ability and the Fermi velocity changes are more sensitive to the changes of lattice constant and buckling angle. We confirm the validation of our tight-binding model by direct first-principles calculations of model systems including hydrogenated monolayer Bi and Sb honeycomb lattices. By offering the generally more intuitive tight binding model, our work presents a new insight into and thus a more in-depth understanding of p_x,y-orbital Dirac states in honeycomb lattices. This will be useful in guiding the Dirac states engineering across the materials and in pursuit of the applications of this family of materials in nanoelectronics.